Spin unrestricted linear scaling electronic structure theory 
and its application to magnetic carbon doped BN nanotubes 
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We present an extension of density matrix based linear scaling electronic structure theory to 
incorporate spin degrees of freedom. When the spin multiplicity of the system can be predetermined, 
the generalization of the existing linear scaling methods to spin unrestricted cases is straightforward. 
However, without calculations it is hard to determine the spin multiplicity of some complex systems, 
such as, many magnetic nanostuctures, some inorganic or bioinorganic molecules. Here we give a 
general prescription to obtain the spin-unrestricted ground state of open shell systems. Our methods 
are implemented into the linear scaling trace-correcting density matrix purification algorithm. The 
numerical atomic orbital basis, rather than the commonly adopted Gaussian basis functions is used. 
The test systems include O2 molecule, and magnetic carbon doped BN(5,5) and BN(7,6) nanotubes. 
Using the newly developed method, we find the magnetic moments in carbon doped BN nanotubes 
couple antiferromagnetically with each other. Our results suggest that the linear scaling spin- 
unrestricted trace-correcting purification method is very powerful to treat large magnetic systems. 

PACS numbers: 



I. INTRODUCTION 

Recently, the spin of the electron has caused a resur- 
gence of interest because it promises a wide variety of 
new devices that combine logic, storage and sensor ap- 
plications. Moreover, these "spintronic" devices might 
lead to quantum computers and quantum communication 
based on electronic solid-state devices. 1 To explain the 
experimental findings and predict novel magnetic proper- 
ties of nanostructures, ah initio electronic structure cal- 
culations on these magnetic systems are indispensable. 
However, ah initio electronic structure calculations are 
usually limited to small and medium size molecular sys- 
tems. The obstacle lies in rapidly increasing computa- 
tional costs as the systems become larger and more com- 
plex. Usually, the molecules which possess novel mag- 
netic properties include few hundreds or thousands of 
atoms. They are too large to be calculated by conven- 
tional ah initio or density functional theory. For instance, 
single- molecule magnets (SMMs), - ' which have attracted 
much interest for the quantum tunneling of magnetiza- 
tion in such systems, are typically very large molecules, 
e.g., [Mn 12 Oi 2 (0 2 CC 6 F 5 ) 1 6(H 2 0) 4 ]-, :i a SMM of Mn 12 
family, is composed by 260 atoms. Extensively studies 
on such systems clearly necessitate new methods with 
desired computational complexity. 

Linear scaling (O(N)) electronic structure theory in 
combination with tight-binding, self-consistent Hartree- 
Fock (HF) or density functional theory (DFT) has be- 
come a very powerful tool to investigate very large com- 
plex systems, ' such as, the silicon systems, ( ' the giant 
fullerenes, the low-dimension nanomaterials,' and the bi- 
ological molecules. 8 They have been successfully applied 
to calculate the molecular energies, obtain the optimized 



molecular geometries and evaluate the static and dynami- 
cal molecular properties of very large systems. 1 Many lin- 
ear scaling methods have already been developed to build 
the effective Hamiltonian" and to avoid cubic scaling 
Roothaan step by replacing it with diagonalization-frcc 
alternatives. 10 ' 11 ' 12 ' 13 - 14 ' 15 ' 16 ' 17 ' 18 ' 19 ' 20 ' 21 ' 22 ' 23 ' 24 - 25 After 
the algorithms whose computational cost scales asymp- 
totically only linearly with the size of the system 
are available for building the effective Hamiltonian, 
the Roothaan step which updates the occupied spaces 
becomes the rate-determining step in large enough 
self-consistent-field (SCF) calculations. Typically this 
may occur for calculations involving several thou- 
sand basis functions. Two of the main alterna- 
tives may be briefly summarized as follows. One 
may attempt to update the one-particle density ma- 
trix itself, 1 ^' 14 ' 15 ' 16 ' 17 - 18 ' 19 ' 20 ' 21 . 22 ' 23 ' 24 ' 25 rather than the 
molecular orbitals. Recent examples include the im- 
proved Fermi Operator expansion method 20 ' 21 ' 22 and the 
purification projection schemes which was first proposed 
by Palser and Manolopoulos 10 and later was improved 
by Niklasson and his coworkers. 23 ' 24,25 Second, one may 
attempt to obtain localized molecular orbitals, 10,11 ' 12 
rather than the delocalized orbitals that diagonalize the 
Hamiltonian matrix. Some systematic comparisons of 
different approaches have been reported in the context of 
semiempirical electronic structure methods. _(1 

However, to the best of our knowledge, linear scaling 
methods are all exclusively applied to closed shell sys- 
tems. As the first attempt of our efforts in the field, 
we provide a first survey on the possibility of applying 
the linear scaling methods to open shell systems at self- 
consistent electronic structure level. We notice that if 
the spin multiplicity of the insulating open shell systems 
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can be predetermined, almost all linear scaling meth- 
ods are readily applicable. But sometimes, without cal- 
culations it is hard to determine the spin multiplicity 
of some complex systems, such as, some inorganic or 
bioinorganic molecules with the transition metal atoms 
involved. Usually one needs to calculate few states as- 
suming different spin multiplicities. By comparing their 
energies one finally determine the spin degree of the sys- 
tems. It is inconvenient and time-consuming. It therefore 
extremely necessary to develop a method which can au- 
tomatically determine spin multiplicity of large complex 
systems while has low computational complexity. 

In this work, we develop a general method to deter- 
mine automatically the spin multiplicity in the calcula- 
tion and extend the Niklasson's trace-correcting density 
matrix purification algorithm (TC2) to deal with spin un- 
restricted systems. Here, we select the TC2 purification 
algorithm to demonstrate and implement our new spin 
unrestricted linear scaling methods since the TC2 purifi- 
cation is very simple and efficient at both low (< 10%) 
and high (> 90%) occupancies. Interestingly, an intrigu- 
ing density matrix perturbation theory 2 ' based on the 
TC2 purification is proposed and succefully applied to 
calculate electric polarizability recently." 8 Our methods 
are implemented in a DFT program employing localized 
numerical atomic orbitals as basis sets. 29 ' 30 This paper 
is organized as follows: In the next section, we shall 
describe the spin unrestricted linear scaling electronic 
structure theory We illustrate our method by present- 
ing a spin unrestricted version of the TC2 algorithm. 
In Sec. Ill, we describe the details of the implemen- 
tation and perform some test calculations to illustrate 
the Tightness, robustness, and linear scaling behavior of 
our methods. Magnetic carbon doped BN nanotubes are 
studied with the spin unrestricted TC2 method. Our re- 
sults clearly indicate that the magnetic moments in car- 
bon doped BN nanotubes couple antiferromagnetically 
with each other. We discuss the possibility of applying 
the spin unrestricted linear scaling methods to magnetic 
metallic systems in Sec. IV. Finally, our concluding re- 
marks are given in Sec. V. 



II. THEORY 

The most general one-particle reduced density matrix 
may be written as 



p(x, x') = p aa (r, r') + p a p(r, r')+ 
P/3a(r,r') + ppp(r,r'), 



(1) 



here x (x') is a combination of a space coordinate r (r') 
and a spin coordinate s (s'). It is shown by McWeeny 51 
that in any state where the z component of total spin 
is definite, the two components p a p(r,r') and pp a (r,r') 
must vanish, and the first-order reduced density matrix 
has the form: 



for simplicity, we write p aa (r,r') as p a (r,r'), and define 
the similar abbreviation for ppp(r,r'). In case where the 
number of electrons for both spin components are inte- 
gers, 



p$ 



Pa 
PP 

N e = N a + Np, 



Tr(p a ) 
Tr(pp) 



Np 



(3) 



here N a and Np are the number of electrons for spin 
up and spin down components respectively, and N e rep- 
resents the total number of electrons. It indicates that 
if N a and Np arc known prior, one can deal with open 
shell systems without any problem using the linear scal- 
ing methods applied to close shell systems. For example, 
given N a and Np, we can deal with open shell systems 
using LNVD density matrix minimization (DMM), > 15 
simplified density matrix minimization (SDMM), 18 or 
spectral projection approach. As an example, here we 
give a spin unrestricted TC2 algorithm with predeter- 
mined spin multiplicity (PSUTC2) to show how it works: 



(Pa.n) 



Tr(p atn ) > N a 

2p a ,n - Pa,n> Tr{p a , n ) < N a 



Pa,n ' 



(4) 



with p afi = {e N {H a )I - H a )/(e afi - e {H a )) (H a is the 
majority part of the Hamiltonian matrix) where the con- 
stants eo(H a ) and £]\[(H a ) are the lowest and highest 
eigenvalues of H a , respectively. The purification algo- 
rithm for pp is similar to the above step just with a 
replaced by j3. We note that McWeeny has extended 
his density matrix method to deal with open shell sys- 
tems with given orbital occupations,' however not in the 
framework of linear scaling methods. However, some- 
times the above algorithm is inconvenient in practical 
applications since it needs the predefined spin multiplic- 
ity. One shortcoming is that to find the electronic ground 
states, one must carry out several calculations with dif- 
ferent spin multiplicity. Another serious problem is that 
one may encounter fractional occupation in practical cal- 
culations if the spin multiplicity is not the right one with 
no fractional occupation. Under this circumstance, the 
algorithm might not be linear scaling using real space 
representations, or even fail since in this case the density 
matrix is not idempotent. 

Here we give a prescription to solve this problem. We 
define new H and p as: 



H = 



H a 
Hp 



P = 



Pa 

pp 



(5) 



p(x,x') = p a (r,r') + pp(r,r'), 



(2) 



Here the new operators satisfy Hp = pH, pp = p. and 
Tr(p) = N e . Through the newly defined Hamiltonian H, 
the spin multiplicity for the ground state can be deter- 
mined automatically according to the Aufbau' principle. 
It is easy to see that if we find p which minimizes Tr(Hp) 
within the above constraints, p a and pp can be extracted 
easily from p. In fact, the minimization problem can be 
solved by many linear scaling methods for closed shell 
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systems. One difference is that we now deaf with ma- 
trices with the dimension 2Nb, here iVf, is the number of 
basis sets. However, due to the block form of H and p, 
the 2N{, dimension problem can be decomposed into two 
Nb dimension problems. Here we detail this method by 
presenting an algorithm based on TC2 projection algo- 
rithm. The spin unrestricted TC2 (SUTC2) projection 
algorithm is given by this pseudocode: 

subroutine SUTC2(H a , Hp, p a , pp, N e , ErrorLimit) 
estimate e (H a ), e N (H a ),e (Hp), e N (Hp) 
en = max(eN(H a ), £ N (Hp)) 
e = min(e (H a ), s {Hp)) 
Pa.o = (snI — H a )/(eN — £o) 
Pl3,o = (snI - Hp)/ (en - e ) 
while Error > ErrorLimit 
if Tr[p }-N e <0 

Pa,n+1 = 2p a ,n — p an 
Pp,n+1 = 2pp,n - Pp tn 

else 

_ 2 
Pa,n+1 Pa.n 

PP,n+l = P 2 p. n 

end 

estimate Error 
end 

Pa Pa,n 
PP = P0,n ■ 

(6) 

The scheme can be described as follows: First using the 
same scaling factors we normalize H a and Hp to get ini- 
tial matrices p Qj o an d pp.o with all theirs eigenvalues in 
the range [0,1]. Then p a . n and pp, n are updated in the 
same way depending on the sum of the traces. In this 
way, p a and pp is obtained using the same purification 
polynomial. The monotonicity of the purification poly- 
nomial leads to the correct occupations according to the 
'Aufbau' principle. 

Although we illustrate our method by presenting only 
the SUTC2 method, however, our method is quite gen- 
eral and can be easily generalized to many other density 
matrix or localized orbitals based linear scaling meth- 
ods. For instance, the KMG functional 11 can be eas- 
ily adapted to include spin degrees of freedom given the 
chemical potential of the magnetic insulating systems. 



ing localized, Wannier-like orbitals employing the KMG 
functional | J in SIESTA. Unfortunately, the convergence 
of the conjugate gradient (CG) minimization of the elec- 
tronic energy in the first SCF step might be extremely 
slow (up to 2000 CG iterations, compared to 20 in fur- 
ther SCF steps). Another inconvenience is the chemi- 
cal potential must be given prior to conserve the total 
charge, which might be notably difficult for small gap 
systems. ' Thus we implement the robust density ma- 
trix purification methods in SIESTA. Saravanan et al. 
showed that the multiatom blocked sparse matrix multi- 
plications can be much faster than a standard element- 
by-element sparse matrix package and also more effi- 
cient than the atom-blocked sparse matrix algebra. 18 The 
blocking scheme benefits from the use of highly-optimized 
level-3 basic linear algebra subroutines (BLAS) for large 
submatrix multiplications. So in this work, we em- 
ploy the blocked compressed sparse row (BCSR) 18 ' 32 ' 33 
data structure with multiatom blocks for sparse ma- 
trix computations. We use a multiatom blocked sparse 
matrix multiplications with dropping (filtering) of mul- 
tiatom blocks with the Frobcnious norm below a numeri- 
cal thrcshold(10~ 4 — 10~ 6 ) to obtain energy as accurate as 
1 mHatrcc. Relative to the cutoff approach, a major ad- 
vantage of threshold metered sparse linear algebra is that 
it avoids discontinuities in the potential energy surface 
associated with atoms moving in and out of the cutoff ra- 
dius. We work in an orthogonal representation though in 
siesta H is evaluated in the nonorthogonal basis of atomic 
orbitals (AO). We achieve this by transforming the AO 
Hamiltonian matrix Hao to an orthonormal basis using 
H = ZHaoZ T an d obtaining the AO density matrix pAO 
using pao = Z T pZ, where the inverse factor Z = L^ 1 , 
and L is the Cholesky factor for which S = LL T . The 
inverse factor Z is obtained directly using the state of the 
art blocked approximate inverse (AINV) algorithm. 3 ,35 
As for the force and stress calculations, only the orthog- 
onality parts one must care about in our density matrix 
implementation. The orthogonality force and stress re- 
quire the energy-density matrix E = paoHaoS^ 1 ■ In 
our implementation, the energy-density matrix is calcu- 
lated as E = {{paoHao)Z t )Z and it is only calculated 
when the SCF reaches its convergence. 



III. IMPLEMENTATION AND RESULTS 

A. Implementation 

Both schemes developed in this work for dealing with 
open shell systems are implemented in SIESTA, a stan- 
dard Kohn-Sham density functional program using norm- 
conserving pseudopotentials and numerical atomic or- 
bitals as basis sets. 21 ' In SIESTA, periodic boundary con- 
ditions are employed to simulate both isolated and peri- 
odic systems. Here since we aim at large systems, r-point 
sampling is used. There is a linear scaling solver us- 



All calculations reported here employ the local density 
approximation (LDA). i(l And no structural optimizations 
are performed in all following calculations to save com- 
puting resource. First we validate our methods and our 
implementation by calculate oxygen molecule with the 
well known triplet ground state. As expected, by speci- 
fying the spin triplet state, using the PSUTC2 method, 
we can get the same energy for the triplet state as that 
from the diagonalization calculation. And the SUTC2 
method can also give the triplet ground state with the 
same energy for oxygen molecule without given the elec- 
tronic occupation. 
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B. Performance of the method 

A recent DFT calculation indicates that carbon doping 
induces spontaneous magnetization in BN nanotubes. 
The density of states (DOS) for both carbon substituted 
BN(5,5) nanotube and carbon substituted BN(9,0) nan- 
otube display the insulating or semiconducting behavior. 
So carbon doping BN nanotubes serve as ideal systems 
to test our methods. We choose BN(5,5) nanotubes with 
a boron atom substituted by a carbon for this purpose. 
Both PSUTC2 and SUTC2 methods are tested. In the 
PSUTC2 method, the magnetic moment of these sys- 
tems is fixed to 1 hb- For these systems, we find that 
the PSUTC2 method is faster by 25% than the SUTC2 
method. In TC2 methods, the efficiency is determined 
by the magnitude of the band gap. 23 In the PSUTC2 
method, the convergence for the spin up density matrix is 
determined by the band gap of the majority part, and so 
does for the spin down component. However, the conver- 
gence in the SUTC2 method depends on the magnitude 
of the system band gap, which is always smaller in fer- 
romagnetic (FM) systems. For antiferromagnetic (AFM) 
systems, the system band gap are the same as the band 
gaps for both spin components, and the SUTC2 method 
should be as efficient as the PSUTC2 method. So, the 
PSUTC2 method is prefered for systems where the spin 
multiplicity is known prior. For systems without knowl- 
edge for the spin multiplicity, one should use the SUTC2 
method to account for the magnetic moment of the sys- 
tems. Furthermore, if we assume the spin multiplicity is 
not changed in all SCF cycles, we can combine the effi- 
ciency of the PSUTC2 method and the robustness of the 
SUTC2 method by using SUTC2 method in the first SCF 
cycle and using PSUTC2 method in the following SCF 
cycles. Fig. 1 shows the CPU time per SCF cycle for 
the diagonalization and PSUTC2 methods. Three dif- 
ferent basis sets (single-^ (SZ), double-^ (DZ), double-^ 
plus polarization functions (DZP)) arc employed. The 
diagonalization method clearly shows a 0(N 3 ) scaling, 
in contrast sharply to the linear scaling behavior of the 
PSUTC2 method for all basis sets. The critical system 
size where the PSUTC2 method is faster than the di- 
agonalization method is dependent on the basis set em- 
ployed. For instance, the critical system size is 320, 300, 
and 200 for SZ, DZ and DZP basis sets respectively. 



C. Applications to magnetic carbon doped BN(5,5) 
and BN(7,6) nanotubes 

Following the recent discovery of fcrromagnctism at 
room temperature in an all-carbon system consisting of 
polymerized C60, 38 there has been increased interest in 
magnetism in metal-free systems. Although previous ex- 
periment indicated a preference for zig-zag and near zig- 
zag BN tubes, 39 a very recent high-resolution electron 
diffraction study on BN nanotubes grown in a carbon- 
free chemical vapor deposition process revealed a disper- 



sion of the chiral angles a in the ranges of 10° < a < 15° 
and 25° < a < 30°. 40 Chiral BN nanotubes usually con- 
tain large number of atoms in a unit cell, e.g., BN(7,6) 
nanotube have 508 atoms in a unit cell with chiral angle 
a = 27.46°. These nanotubes are difficult to be treated 
using traditional methods. Here we conduct a linear scal- 
ing spin unrestricted calculation on BN(7,6) nanotube 
with a boron atom substituted by a carbon. The DZ basis 
set is employed for this purpose. For better comparison 
with carbon doped BN(5,5) nanotube, we show the struc- 
ture and calculated spin density for both carbon doped 
BN(5,5) nanotube with 500 atoms per cell and and car- 
bon doped BN(7,6) nanotube. Clearly, the distribution 
of spin density for both carbon doped nanotubes is ba- 
sically the same: The spin density is mainly contributed 
by carbon 2p z orbital and the nitrogen atoms near car- 
bon have some small magnetic moments. To learn more 
knowledge about the electronic and magnetic properties 
of the system, we obtain the density of states (DOS) 
by a non-SCF calculation using diagonalization. Fig. 3 
shows the total DOS (TDOS) as well as carbon partial 
DOS (PDOS) for both carbon doped BN(5,5) nanotube 
with 500 atoms per cell and carbon doped BN(7,6) nan- 
otube. We can see that in the plotting energy range, in 
both cases carbon has sizable contribution only to the two 
states around the Fermi level. The shape of the DOS is 
very similar for both carbon doped BN nanotubes. Also 
these DOS resemble that of carbon doped BN(5,5) nan- 
otube with 80 atoms per cell provided by Wu et al. 
except that the band gap in our case is smaller possibly 
due to the LDA functional we employed. 

Previous calculations indicated there exist magnetic 
moments in carbon doped BN nanotubes and our cal- 
culations confirm it. Then, how the local magnetic mo- 
ments in carbon doped BN nanotubes couple with each 
other: ferromagnetically or antiferromagnetically? Here 
we conduct some calculations on BN(7,6) nanotube with 
two boron atoms substituted by two carbon atoms to ad- 
dress this issue. Six different configurations (A, B, C, D, 
E, and F) are considered, as shown in Fig. 4(a). The total 
energy of both AFM and FM states for these six carbon 
doped BN nanotubes is shown in Fig. 4(b). Generally, the 
relatively energy difference converges much faster along 
with the cutoff radius for the density matrix than the 
absolute total energy. In fact, some test calculations in- 
dicate that the deviation of the energy difference from 
the exact value is smaller than 10 meV. As shown in 
Fig. 4(b), structure B with AFM state is the most sta- 
ble configuration among these six carbon doped BN(7,6) 
nanotubes. Moreover, we can see that when two carbon 
atoms are in the same hexagonal ring, the AFM state is 
more favorable over the FM state. In addition, the en- 
ergy difference for B is larger than that for A possibly 
due to the fact that the two carbon 2p z orbitals in B are 
more parallel and thus the superexchange AFM interac- 
tion in B is larger. If the distance between two carbon 
atoms is larger than 3 A, i.e., two carbon atoms are not 
in the same hexagonal ring, the coupling between two 
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magnetic moments is negligibly small. So from Fig. 4(b), 
we can conclude that the magnetic moments in carbon 
doped BN nanotubes couple antifcrromagnctically with 
each other. 



IV. DISCUSSION 

Recently, developing linear scaling methods to deal 
with metallic systems has attracted much interest. Baer 
and Head-Gordon developed an energy rcnormalization- 
group method, with which the computational effort scales 
near linearly with system size even when the density ma- 
trix is highly nonlocal. 41 Goedecker et al. showed that in 
addition to the real-space localization of the density ma- 
trix, there is also a localization in Fourier space. Within 
the multircsolution wavelet representation it was shown 
how the sparsity of the density matrix is preserved for lo- 
calized insulating systems as well as for itinerant metallic 
systems. Combining such techniques with our methods, 
O(N) calculations on metallic magnetic systems might 
be possible since the number of matrix multiplications in 
the TC2 purification, essentially is independent of sys- 
tem size even for metallic systems. 2 ' 5 Although, there 
are no unfilled shells in spin unrestricted HF theory, as 
proved by Bach et al., i3 however, for spin unrestricted 
DFT methods, sometimes fractional occupation is possi- 
ble. In principle, the fractional occupation problem can 
be treated using PM 19 or TRS4 J ! purification. So the 
generalized spin unrestricted form of PM or TRS4 purifi- 
cation can treat these systems. 



V. CONCLUSIONS 

To conclude, we give a first survey on applying linear 
scaling electronic structure methods to spin unrestricted 
systems. Two methods are proposed to deal with systems 
with or without predetermined spin multiplicity respec- 
tively. We demonstrate our methods by detailing the 
PSUTC2 and SUTC2 projection algorithms. The cur- 
rent methods have been implemented in a Kohn-Sham 
density functional program employing numerical atomic 
orbitals as basis sets. We apply our method to deal with 
magnetic carbon doped BN nanotubes. Carbon doped 
BN(7,6) nantoubc has similar magnetic properties as car- 
bon doped BN(5,5) nantoube. Moreover, FM coupling is 
unfavorable for carbon doped BN nantoubes. The results 
suggest that our methods pave the way for carrying out 
linear scaling calculations on spin unrestricted systems, 
such as magnetic nanostructures. 
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FIG. 1: CPU time per SCF cycle for carbon doped BN(5,5) 
nanotube using the PSUTC2 method and traditional diago- 
nalization method with different basis sets. All calculations 
were carried out on a 1.5 GHz Itanium 2 CPU processor run- 
ning RedHat Linux Advanced Server V2.1. 
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FIG. 2: Structures for (a) carbon doped BN(5,5) nanotube 
with 500 atoms per cell and (b) carbon doped BN(7,6) nan- 
otube. The regions between the vertical lines indicate the 
unit cells. Spin densities of the two nanotubes are shown in 
(c) and (d) respectively. 
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FIG. 3: Majority and minority spin DOS of (a) carbon 
doped BN(5,5) nanotube with 500 atoms per cell and (b) 
carbon doped BN(7,6) nanotube. TDOS and carbon PDOS 
are shown in black and red solid lines respectively. The Fermi 
energy is indicated by the vertical dashed line. 
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FIG. 4: (a) Six different structures A, B, C, D, E, and F 
for BN(7,6) tubes with two boron atoms substituted by two 
carbon atoms, (b) shows the total energy for both AFM and 
FM states for different carbon doped BN(7,6) tubes. The 
horizontal axis is the nearest distance between two carbon 
atoms. 



